3. Linear models

---

• 필요 패키지 : car, ISLR, MASS, corrplot
• 필요 데이터 : Boston{MASS}, Carseats{ISLR}
• 주요 함수 : lm{stats}, vif{car}

---

3.1 상관분석 correlation analysis

• 연속형인 두 변수 사이의 종속성을 측정하는 도구인 상관계수(correlation coefficient)를 기초로 해 시행하는 분석 절차를 상관분석이라 함

• 두 변수 사이의 함수 관계를 규명하는 회귀분석의 전 단계에 산점도 분석과 함께 반드시 시행되는 절차임

• 상관계수는 Karl Pearson (1857 – 1936)이 제안한 개념으로 정의는 다음과 같음

\[ \rho = \frac{\textsf{Cov}(X, Y)}{\textsf{sd}(X)\textsf{sd}(Y)} \]

• 이 측도는 두 연속형 변수 사이의 선형 종속성(linear dependency)의 정도를 측정

• \(-1\) ~ \(+1\) 사이의 값을 가지며, 0 근처이면 선형종속성이 유의하지 않음을, +1 또는 -1이면 완벽한 선형종속성을 가짐을 의미

Beyond 선형종속성

• 스피어만(Spearman)의 순위상관계수는 변수의 순위 사이의 상관계수를 계산한 것으로 선형 종속성을 넘어 단조적(monotone) 종속성을 측정해 줌

• 켄달(Kendall)의 타우(\(\tau\))는 데이터에서 두 자료점을 선택해 concordant/discordant 여부를 따져서 concordant pair와 discordant pair 개수의 차이를 이용해 변수 간 연관성을 측정하는 도구임

round(cor(x1, x7, method = "spearman"), 4)

## [1] 0.7189

round(cor(x1, x7, method = "kendall"), 4)

## [1] 0.5167

round(cor(x10, x11, method = "spearman"), 4)

## [1] 0.5217

round(cor(x10, x11, method = "kendall"), 4)

## [1] 0.3453

예제 데이터 분석

• iris 데이터

pairs(iris[,1:4], 
      main = "Anderson's Iris Data -- 3 species", pch = 21, 
      bg = c("red", "green3", "blue")[unclass(iris$Species)])

• 이 데이터에서 setosa 품종 데이터만 추출해 변수들 간의 상관성을 분석

Setosa <- subset(iris, Species == "setosa", select = -Species)
round(cor(Setosa), 4)

##              Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length       1.0000      0.7425       0.2672      0.2781
## Sepal.Width        0.7425      1.0000       0.1777      0.2328
## Petal.Length       0.2672      0.1777       1.0000      0.3316
## Petal.Width        0.2781      0.2328       0.3316      1.0000

• 스피어만의 순위상관계수와 켄달의 타우

round(cor(Setosa, method = "spearman"), 4)

##              Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length       1.0000      0.7553       0.2789      0.2995
## Sepal.Width        0.7553      1.0000       0.1799      0.2865
## Petal.Length       0.2789      0.1799       1.0000      0.2711
## Petal.Width        0.2995      0.2865       0.2711      1.0000

round(cor(Setosa, method = "kendall"), 4)

##              Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length       1.0000      0.5973       0.2173      0.2311
## Sepal.Width        0.5973      1.0000       0.1426      0.2343
## Petal.Length       0.2173      0.1426       1.0000      0.2217
## Petal.Width        0.2311      0.2343       0.2217      1.0000

• 두 변수 사이의 상관관계가 있는 지 여부를 검정할 수도 있음.

  • 이때 귀무가설은 \(H_0: \rho = 0\)

with(Setosa, cor.test(Sepal.Width, Petal.Width))

## 
##  Pearson's product-moment correlation
## 
## data:  Sepal.Width and Petal.Width
## t = 1.6581, df = 48, p-value = 0.1038
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.0487543  0.4800023
## sample estimates:
##      cor 
## 0.232752

with(Setosa, cor.test(Sepal.Width, Petal.Width, 
                      method = "spearman", exact = FALSE))

## 
##  Spearman's rank correlation rho
## 
## data:  Sepal.Width and Petal.Width
## S = 14858, p-value = 0.04366
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2865359

with(Setosa, cor.test(Sepal.Width, Petal.Width, 
                      method = "kendall", exact = FALSE))

## 
##  Kendall's rank correlation tau
## 
## data:  Sepal.Width and Petal.Width
## z = 2.0593, p-value = 0.03946
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.2342674

---

3.2 선형회귀모형 Linear regression model

• 회귀분석이란?

  • 변수들 간의 함수적 관계를 밝히려는 통계적 분석 방법으로 영국의 유전학자이자 우생학의 창시자인 갈튼이 1885년에 발표한 논문에서 유래된 지도학습 방법

• 종속변수(dependent variable) 혹은 반응변수(response variable)

  • 다른 변수에 의해 설명되는 변수

  • 보통 측정오차를 수반함

• 독립변수(independent variable) 혹은 설명변수(explanatory variable)

  • 다른 변수에 영향을 주거나 그 변수를 설명하는 변수

  • 보통 측정오차가 없다고 가정

• 회귀모형 \[ y = f(x) + \epsilon \]

• 회귀분석의 활용

  • 여러 통계적인 분석기법 가운데 가장 널리 쓰임

  • 공학, 자연과학, 사회과학 및 의약학 등 모든 학문 분야에서 사용됨

  • 분류 문제도 회귀분석 문제로 볼 수 있음

• 사례 : Advertising data

• Questions

  • 각 매체 광고비와 판매량간에 어떤 관계가 있는가?

  • 관계가 있다면 얼마나 강해 보이는가?

  • 광고비와 판매량 간의 관계가 선형인가?

---

3.3 단순선형회귀 simple linear regression

• 반응변수 \(Y\)와 하나의 설명변수 \(X\)간에 선형함수관계를 가정하여 모형을 적합

• 실제 회귀함수가 선형이 아닐지라도 선형회귀는 해석 및 활용이 간편하다는 장점을 가짐

• 단순선형회귀모형

\[Y = \beta_0 + \beta_1 X + \varepsilon\]

• 기본가정 : 모형의 선형성, 오차항의 정규성, 독립성, 등분산성

• 최소제곱추정(least square estimate)

  • Residual Sum of Squares (RSS)를 최소로 만드는 회귀계수의 값을 구함

  \[\mbox{RSS} = (y_1-\hat\beta_0 - \hat\beta_1 x_1)^2+(y_2-\hat\beta_0 - \hat\beta_1 x_2)^2+\cdots+(y_n-\hat\beta_0 - \hat\beta_1 x_n)^2\]

• Measures of Fit: \(R^2\)

  • 위 회귀모형에 따르면 \(Y\)의 변동의 요인은 설명변수 \(X\)와 오차 \(\epsilon\) 두 가지임

  • \(Y\)의 변동을 설명변수에 의한 변동과 오차에 의한 변동으로 분해할 수 있음 (분산분석, ANOVA)

  • \(R^2\)는 \(Y\)의 변동 중 설명변수 \(X\)에 의한 변동의 상대적 크기를 의미 \[ R^2 = 1 - \frac{\text{RSS}} {\sum(Y_i-\bar{Y})^2} \approx 1 - \frac{\text{Ending Variance}}{\text{Starting Variance}}\]

  • \(R^2\)는 항상 0과 1 사이의 값을 가짐

  • \(R^2 = 0\)은 \(X\)를 이용해 설정했던 모형이 아무런 설명력이 없음을 의미하고, \(R^2 = 1\)은 그 모형으로 데이터를 완벽하게 설명한다는 것을 의미

Example: 모의실험 자료

• 자료 생성을 위한 모형식 \[ y = 2 + 3 x + \varepsilon, \quad \varepsilon \sim N(0, 5^2) \]

set.seed(123)
N <- 50
x <- runif(N, min = -3, max = 3)
y <- 2 + 3*x + rnorm(N, sd = 5)
plot(x, y)
abline(2, 3, col = "navy", lwd = 3)

• 모형 적합

fit.lm <- lm(y ~ x)
plot(x, y)
abline(2, 3, col = "navy", lwd = 3)
abline(print(coef(fit.lm)), col = 2, lwd = 3)

## (Intercept)           x 
##    2.242598    3.318121

legend("topleft", c("True", "Fitted"), col = c("navy", "red"), lwd = c(3, 3))

• 분석 결과 정리

summary(fit.lm)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.2789  -2.7893  -0.3284   2.7463  10.9306 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.2426     0.6665   3.365  0.00151 ** 
## x             3.3181     0.3804   8.723 1.82e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.702 on 48 degrees of freedom
## Multiple R-squared:  0.6132, Adjusted R-squared:  0.6051 
## F-statistic: 76.08 on 1 and 48 DF,  p-value: 1.822e-11

• \(p\)-값의 의미

M <- 10000
beta1 <- rep(NA, M)
for ( r in 1:M ) {
  y <- 2 + 0*x + rnorm(N, sd = 5)
  beta1[r] <- coef(lm(y ~ x))[2]
}

hist(beta1, probability= TRUE, main = "", xlim = c(-3, 5), 
     col = "lightgray", border = "white", nclass = 100,
     xlab = expression(beta[1]))
abline(v = coef(fit.lm)[2], col = 2, lty = 2)
points(coef(fit.lm)[2], 0, pch = 25, bg = "yellow", col = 2)

mean(beta1 >= coef(fit.lm)[2])

## [1] 0

Example: Advertising data

• 데이터

adv <- read.csv(file = "./dat/Advertising.csv")
head(adv)

##   X    TV Radio Newspaper Sales
## 1 1 230.1  37.8      69.2  22.1
## 2 2  44.5  39.3      45.1  10.4
## 3 3  17.2  45.9      69.3   9.3
## 4 4 151.5  41.3      58.5  18.5
## 5 5 180.8  10.8      58.4  12.9
## 6 6   8.7  48.9      75.0   7.2

summary(adv)

##        X                TV             Radio          Newspaper     
##  Min.   :  1.00   Min.   :  0.70   Min.   : 0.000   Min.   :  0.30  
##  1st Qu.: 50.75   1st Qu.: 74.38   1st Qu.: 9.975   1st Qu.: 12.75  
##  Median :100.50   Median :149.75   Median :22.900   Median : 25.75  
##  Mean   :100.50   Mean   :147.04   Mean   :23.264   Mean   : 30.55  
##  3rd Qu.:150.25   3rd Qu.:218.82   3rd Qu.:36.525   3rd Qu.: 45.10  
##  Max.   :200.00   Max.   :296.40   Max.   :49.600   Max.   :114.00  
##      Sales      
##  Min.   : 1.60  
##  1st Qu.:10.38  
##  Median :12.90  
##  Mean   :14.02  
##  3rd Qu.:17.40  
##  Max.   :27.00

par(mfrow = c(1, 3))
with(adv, plot(TV, Sales))
with(adv, plot(Radio, Sales))
with(adv, plot(Newspaper, Sales))

• 반응변수: Sales

• 독립변수: TV

• 회귀모형식: \[ \textsf{Sales} = \beta_0 + \beta_1 \cdot \textsf{TV} + \varepsilon \]

• 모형 적합

  • 함수 lm()을 사용하면 간단

fit.lm <- lm(Sales ~ TV, data = adv)
with(adv, plot(TV, Sales))
abline(coef(fit.lm), col = "navy", lwd = 3)

• 분석 결과 정리

  • 위에서 적합된 결과를 범용함수인 summary()를 이용해 요약

summary(fit.lm)

## 
## Call:
## lm(formula = Sales ~ TV, data = adv)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.3860 -1.9545 -0.1913  2.0671  7.2124 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 7.032594   0.457843   15.36   <2e-16 ***
## TV          0.047537   0.002691   17.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.259 on 198 degrees of freedom
## Multiple R-squared:  0.6119, Adjusted R-squared:  0.6099 
## F-statistic: 312.1 on 1 and 198 DF,  p-value: < 2.2e-16

---

3.4 다중회귀분석 multiple linear regression

• 설명변수가 TV 하나인 것보다는 Radio, Newspaper까지 다 사용하는 것이 어떨까?

• 다중회귀모형식 \[Y_i = \beta_{0} + \beta_{1}X_{1} + \beta_{2}X_{2} + \cdots + \beta_{p}X_{p} + \varepsilon\]

  • \(\beta_0\)는 절편

  • \(\beta_j\)는 \(X_j\)가 한 단위 증가하고 나머지 \(X\) 변수의 값은 일정하게 유지될 때, \(Y\)의 평균 증가량에 해당 (partial regression coefficients)

• 최소제곱추정: 다음 식을 최소로 하는 \(\beta\) 계수

\[\text{MSE} = \frac{1}{n} \sum_{i=1}^n (Y_i-\hat{Y})^2 = \frac{1}{n} \sum_{i=1}^n (Y_i-\hat{\beta_{0}}-\hat{\beta_{1}}X_{1}-\cdots-\hat{\beta_{p}}X_{p})^2\]

fit.lm <- lm(Sales ~ TV + Radio + Newspaper, data = adv)

• \(\beta_j\)가 0 인가?… 에 대한 답은 가설검정을 통해서 얻을 수 있음

  • \(\beta_j = 0\) 이면 \(X_j\)를 모형에 포함시킬 이유가 없음

  • 각 계수별 \(t\)-test 결과 확인

summary(fit.lm)

## 
## Call:
## lm(formula = Sales ~ TV + Radio + Newspaper, data = adv)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.8277 -0.8908  0.2418  1.1893  2.8292 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.938889   0.311908   9.422   <2e-16 ***
## TV           0.045765   0.001395  32.809   <2e-16 ***
## Radio        0.188530   0.008611  21.893   <2e-16 ***
## Newspaper   -0.001037   0.005871  -0.177     0.86    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.686 on 196 degrees of freedom
## Multiple R-squared:  0.8972, Adjusted R-squared:  0.8956 
## F-statistic: 570.3 on 3 and 196 DF,  p-value: < 2.2e-16

• \(X\) 변수들 중에서 적어도 하나는 유용한 예측 변수라고 확신할 수 있는가?

  • 이 질문은 결국 \(\beta_1\) = \(\beta_2\) = \(\dots\) = \(\beta_p\) = 0 인가?… 와 같은 것임

  • \(F\)-test 결과 확인

• 모형 수정

  • Newspaper는 유의하지 않은 것으로 나타남

  • Newspaper를 제외한 나머지 두 개의 변수만을 포함한 회귀모형

fit.lm <- lm(Sales ~ TV + Radio, data = adv)
summary(fit.lm)

## 
## Call:
## lm(formula = Sales ~ TV + Radio, data = adv)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.7977 -0.8752  0.2422  1.1708  2.8328 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.92110    0.29449   9.919   <2e-16 ***
## TV           0.04575    0.00139  32.909   <2e-16 ***
## Radio        0.18799    0.00804  23.382   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.681 on 197 degrees of freedom
## Multiple R-squared:  0.8972, Adjusted R-squared:  0.8962 
## F-statistic: 859.6 on 2 and 197 DF,  p-value: < 2.2e-16

• 예측 정확도 평가

  • 예측식: \[\textsf{Sales} = 2.9211 + 0.0458 \cdot \textsf{TV} + 0.1880 \cdot \textsf{Radio}\]

Sales.pred <- predict(fit.lm)
plot(adv$Sales, Sales.pred, xlab = "Actual", ylab = "Predicted")
abline(0, 1, col = 2, lty = 2)

mean((adv$Sales - Sales.pred)^2)  # mean(fit.lm2$residuals^2)

## [1] 2.78457

• TV 광고비로 100, Radio 광고비로 45를 사용한 경우 Sales 값 예측

\[ \textsf{Sales} = 2.9211 + 0.0458 \cdot 100 + 0.1880 \cdot 45 \]

newx <- data.frame(TV = 100, Radio = 45)
predict(fit.lm, newx)

##        1 
## 15.95632

How to choose important variables?

• 확실한 방법은 all possible regression, 즉 모든 가능한 경우를 고려해 최적 모형을 찾는 것임

• 그러나 too many!!! (\(2^p\) with \(p\) predictors)

• Forward selection (전진선택법) :

  • 절편항만 포함하는 Null model로부터 출발하여 RSS 기준으로 변수를 하나씩 추가해 가는 방법

  • 남은 변수가 모두 유의하지 않을 때까지 반복해서 시행

• Backward elimination (후진소거법) :

  • 모든 변수가 들어있는 Full model로부터 출발하여 가장 큰 p-value를 갖는 유의하지 않은 변수를 제거해 가는 방법

  • 남은 변수가 모두 유의할 때까지 반복해서 시행

• Stepwise regression (단계별 선택법) :

  • 변수를 하나씩 추가하면서 새로 들어온 변수로 인해 유의하지 않은 변수는 제거해 가는 방법

  • 새로운 변수의 추가가 없으면 중지

  • 함수 step()을 이용하면 간단

step(fit.lm)

## Start:  AIC=210.82
## Sales ~ TV + Radio
## 
##         Df Sum of Sq    RSS    AIC
## <none>                556.9 210.82
## - Radio  1    1545.6 2102.5 474.52
## - TV     1    3061.6 3618.5 583.10

## 
## Call:
## lm(formula = Sales ~ TV + Radio, data = adv)
## 
## Coefficients:
## (Intercept)           TV        Radio  
##     2.92110      0.04575      0.18799

• 최적의 모델을 선택하기 위한 기준

  • \(C_p\), AIC, BIC, adjusted-\(R^2\), CV 등등

3.5 질적 설명변수 Qualitative Predictors

• 예제 데이터로 ISLR 패키지의 Credit데이터를 사용

library(ISLR)
summary(Credit)

##        ID            Income           Limit           Rating     
##  Min.   :  1.0   Min.   : 10.35   Min.   :  855   Min.   : 93.0  
##  1st Qu.:100.8   1st Qu.: 21.01   1st Qu.: 3088   1st Qu.:247.2  
##  Median :200.5   Median : 33.12   Median : 4622   Median :344.0  
##  Mean   :200.5   Mean   : 45.22   Mean   : 4736   Mean   :354.9  
##  3rd Qu.:300.2   3rd Qu.: 57.47   3rd Qu.: 5873   3rd Qu.:437.2  
##  Max.   :400.0   Max.   :186.63   Max.   :13913   Max.   :982.0  
##      Cards            Age          Education        Gender    Student  
##  Min.   :1.000   Min.   :23.00   Min.   : 5.00    Male :193   No :360  
##  1st Qu.:2.000   1st Qu.:41.75   1st Qu.:11.00   Female:207   Yes: 40  
##  Median :3.000   Median :56.00   Median :14.00                         
##  Mean   :2.958   Mean   :55.67   Mean   :13.45                         
##  3rd Qu.:4.000   3rd Qu.:70.00   3rd Qu.:16.00                         
##  Max.   :9.000   Max.   :98.00   Max.   :20.00                         
##  Married              Ethnicity      Balance       
##  No :155   African American: 99   Min.   :   0.00  
##  Yes:245   Asian           :102   1st Qu.:  68.75  
##            Caucasian       :199   Median : 459.50  
##                                   Mean   : 520.01  
##                                   3rd Qu.: 863.00  
##                                   Max.   :1999.00

• Gender와 같은 범주형변수를 어떻게 회귀식에 포함시킬 것인가?

• Dummy variable (가변수)로 데이터 코딩

  • Gender와 같은 이항변수의 경우 Male이면 0, Female이면 1로 입력

  • Ethnicity 같은 다항변수의 경우에는 one-hot 인코딩: African American는 (1, 0), Asian은 (0, 1)

3.5.1 모형의 해석 Interpretation

• Balance를 반응변수로, Income과 Gender을 설명변수로 사용

• 이항변수인 Gender를 다음과 같이 코딩

\[\textsf{Gender} = \begin{cases} 0 & \quad \text{if male} \\ 1 & \quad \text{if female}\\ \end{cases}\]

• 회귀식은 다음과 같이 표현됨

\[\textsf{Balance} = \beta_{0} + \beta_{1} \cdot \textsf{Income} + \beta_{2} \cdot \textsf{Gender} = \begin{cases} \beta_{0} + \beta_{1} \cdot \textsf{Income} & \quad \text{if male} \\ \beta_{0} + \beta_{1} \cdot \textsf{Income} + \beta_{2} & \quad \text{if female}\\ \end{cases}\]

fit.lm <- lm(Balance ~ Income + Gender, data = Credit)
summary(fit.lm)

## 
## Call:
## lm(formula = Balance ~ Income + Gender, data = Credit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -791.23 -351.34  -51.57  328.18 1112.87 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  233.7663    39.5322   5.913 7.24e-09 ***
## Income         6.0521     0.5799  10.437  < 2e-16 ***
## GenderFemale  24.3108    40.8470   0.595    0.552    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 408.2 on 397 degrees of freedom
## Multiple R-squared:  0.2157, Adjusted R-squared:  0.2117 
## F-statistic: 54.58 on 2 and 397 DF,  p-value: < 2.2e-16

• 여기서, \(\beta_2\)의 추정치는 남성(Male) 대비 여성(Females)의 Balance 초과값의 평균을 의미

3.5.2 다항 변수 More than two levels

• 예를 들어, 3가지 값(Asian, Caucasian, African American)을 갖는 인종 변수의 경우는 다음 2개의 가변수를 사용하고, 추정하는 회귀식은 아래와 같음

\[x_{i1} = \begin{cases} 1 & \quad \text{if Asian} \\ 0 & \quad \text{if not Asian}\\ \end{cases}\]

\[x_{i2} = \begin{cases} 1 & \quad \text{if Caucasian} \\ 0 & \quad \text{if not Caucasian}\\ \end{cases}\]

\[Y_i = \beta_{0} + \beta_{1}\, x_{i1} + \beta_{2}\, x_{i2} + \varepsilon_i = \begin{cases} \beta_{0} + \beta_{1} + \varepsilon_i & \quad \text{if Asian} \\ \beta_{0} + \beta_{2} + \varepsilon_i & \quad \text{if Caucasian}\\ \beta_{0} + \varepsilon_i & \quad \text{if African American}\\ \end{cases}\]

3.6 교호작용 Interaction

3.6.1 양적변수 사이의 교호작용

• \(X_1\)의 증가로 인한 \(Y\)의 (증가 혹은 감소와 같은) 효과가 \(X_2\)값의 수준에 따라 달라질 때, interaction이 존재한다고 함

• 예제: Advertising 데이터

  • TV와 Radio의 광고는 모두 판매량을 증가시킴

  • 두 가지 미디어 광고에 돈을 사용하는 것이 한 가지 광고에만 사용하는 것보다 판매량을 더 많이 증가시킴. Why?

\[\textsf{Sales} = \beta_{0} + \beta_{1} \times \textsf{TV} + \beta_{2} \times \textsf{Radio} + \beta_{12} \times \textsf{TV} \times \textsf{Radio}\]

\[\textsf{Sales} = \beta_{0} + (\beta_{1} + \beta_{12} \times \textsf{Radio}) \times \textsf{TV} + \beta_{2} \times \textsf{Radio}\]

\[\textsf{Sales} = \beta_{0} + (\beta_{2} + \beta_{12} \times \textsf{TV}) \times \textsf{Radio} + \beta_{2} \times \textsf{TV}\]

fit.lm <- lm(Sales ~ TV * Radio, data = adv)
summary(fit.lm)

## 
## Call:
## lm(formula = Sales ~ TV * Radio, data = adv)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3366 -0.4028  0.1831  0.5948  1.5246 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.750e+00  2.479e-01  27.233   <2e-16 ***
## TV          1.910e-02  1.504e-03  12.699   <2e-16 ***
## Radio       2.886e-02  8.905e-03   3.241   0.0014 ** 
## TV:Radio    1.086e-03  5.242e-05  20.727   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9435 on 196 degrees of freedom
## Multiple R-squared:  0.9678, Adjusted R-squared:  0.9673 
## F-statistic:  1963 on 3 and 196 DF,  p-value: < 2.2e-16

• TV에 돈을 더 사용하는 것은 평균적인 판매량을 0.0191 + 0.0011 \(\times\) Radio만큼의 빠르기로 증가시킴

• Radio에 돈을 더 사용하는 것은 평균적인 판매량을 0.0289 + 0.0011 \(\times\) TV으로 증가시킴

library(plot3D)

x <- adv$TV
y <- adv$Radio
z <- adv$Sales
fit.lm1 <- lm(z ~ x + y)
fit.lm2 <- lm(z ~ x * y)

grid.lines <- 200
x.pred <- seq(min(x), max(x), length.out = grid.lines)
y.pred <- seq(min(y), max(y), length.out = grid.lines)
xy <- expand.grid(x = x.pred, y = y.pred)

z1.pred <- matrix(predict(fit.lm1, newdata = xy), 
                  nrow = grid.lines, ncol = grid.lines)
fitpoints1 <- predict(fit.lm1)

z2.pred <- matrix(predict(fit.lm2, newdata = xy), 
                  nrow = grid.lines, ncol = grid.lines)
fitpoints2 <- predict(fit.lm2)


scatter3D(x, y, z, pch = 18, cex = 2,  
          theta = 50, phi = 16, ticktype = "detailed",
          xlab = "TV", ylab = "Radio", zlab = "Sales",  
          surf = list(x = x.pred, y = y.pred, z = z1.pred, 
                      fit = fitpoints1, col = "white", shade = 0.05), 
          main = "Advertising: w/o interaction", col = "red", alpha = 0.5)

scatter3D(x, y, z, pch = 18, cex = 2,  
          theta = 50, phi = 16, ticktype = "detailed",
          xlab = "TV", ylab = "Radio", zlab = "Sales",  
          surf = list(x = x.pred, y = y.pred, z = z2.pred, 
                      fit = fitpoints2, col = "white", shade = 0.05), 
          main = "Advertising: with interaction", col = "red", alpha = 0.5)

---

• Note : Hierarchical principle states that if we include an interaction in a model, we should also include the main effects, even if the p-values associated with their coefficients are not significant.

3.6.2 질적변수와의 교호작용

• 예제: Credit 데이터 \[ \textsf{Balance} = \beta_{0} + \beta_{1} \times \textsf{Income} + \beta_{2} \times \textsf{Student} \] vs. \[ \textsf{Balance} = \beta_{0} + \beta_{1} \times \textsf{Income} + \beta_{2} \times \textsf{Student} + \beta_{12}\times \textsf{Income}\times\textsf{Student} = \begin{cases} \beta_{0} + \beta_{1} \times \textsf{Income} & \quad \text{if}\,\textsf{Non-Student} \\ (\beta_{0} + \beta_{2}) + (\beta_{1} + \beta_{12})\times\textsf{Income} & \quad \text{if}\, \textsf{Student}\\ \end{cases} \]

fit.lm1 <- lm(Balance ~ Income + Student, data = Credit)
fit.lm2 <- lm(Balance ~ Income * Student, data = Credit)
summary(fit.lm1)

## 
## Call:
## lm(formula = Balance ~ Income + Student, data = Credit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -762.37 -331.38  -45.04  323.60  818.28 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 211.1430    32.4572   6.505 2.34e-10 ***
## Income        5.9843     0.5566  10.751  < 2e-16 ***
## StudentYes  382.6705    65.3108   5.859 9.78e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 391.8 on 397 degrees of freedom
## Multiple R-squared:  0.2775, Adjusted R-squared:  0.2738 
## F-statistic: 76.22 on 2 and 397 DF,  p-value: < 2.2e-16

summary(fit.lm2)

## 
## Call:
## lm(formula = Balance ~ Income * Student, data = Credit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -773.39 -325.70  -41.13  321.65  814.04 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       200.6232    33.6984   5.953 5.79e-09 ***
## Income              6.2182     0.5921  10.502  < 2e-16 ***
## StudentYes        476.6758   104.3512   4.568 6.59e-06 ***
## Income:StudentYes  -1.9992     1.7313  -1.155    0.249    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 391.6 on 396 degrees of freedom
## Multiple R-squared:  0.2799, Adjusted R-squared:  0.2744 
## F-statistic:  51.3 on 3 and 396 DF,  p-value: < 2.2e-16

with(Credit, plot(Income, Balance, type = "n", main = "Model without interaction terms"))
beta <- coef(fit.lm1)
abline(beta[1], beta[2], lwd = 3)
abline(beta[1] + beta[3], beta[2], lwd = 3, col = 2)
legend("topleft", c("Student", "Non-student"), col = c(2, 1), lwd = c(3, 3))

with(Credit, plot(Income, Balance, type = "n", main = "Model with interaction terms"))
beta <- coef(fit.lm2)
abline(beta[1], beta[2], lwd = 3)
abline(beta[1] + beta[3], beta[2] + beta[4], lwd = 3, col = 2)
legend("topleft", c("Student", "Non-student"), col = c(2, 1), lwd = c(3, 3))

---

3.7 Potential Fit Problems

3.7.1 Non-linearity of the data

• 다항회귀 polynomial regression

  • 회귀식이 설명변수의 일차 선형식이 아닌 경우

with(Auto, plot(horsepower, mpg, col = "gray", main = "Auto dataset"))

fit.lm1 <- lm(mpg ~ horsepower, data = Auto)
fit.lm2 <- lm(mpg ~ horsepower + I(horsepower^2), data = Auto)

x <- seq(from = min(Auto$horsepower), to = max(Auto$horsepower), by = 0.0025)
with(Auto, plot(horsepower, mpg, col = "gray"))
pred1 <- predict(fit.lm1, newdata = data.frame(horsepower = x))
lines(x, pred1, col = "orange", lwd = 3)
pred2 <- predict(fit.lm2, newdata = data.frame(horsepower = x, x^2))
lines(x, pred2, col = "darkgreen", lwd = 3)
legend("topright", c("Linear", "Quadratic"), col = c("orange", "darkgreen"), lwd = c(3, 3))

• Note: 비선형회귀 nonlinear regression

  • 회귀식이 설명변수의 비선형함수인 경우

3.7.2 Non-constant variance of error terms

• 최소제곱추정법은 오차항의 등분산성을 가정할 때 최고의 방법임

• 그러나 종종 이 가정이 깨지는 경우가 발생

fit.lm <- lm(mpg ~ horsepower + I(horsepower^2), data = Auto)
summary(fit.lm)

## 
## Call:
## lm(formula = mpg ~ horsepower + I(horsepower^2), data = Auto)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.7135  -2.5943  -0.0859   2.2868  15.8961 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     56.9000997  1.8004268   31.60   <2e-16 ***
## horsepower      -0.4661896  0.0311246  -14.98   <2e-16 ***
## I(horsepower^2)  0.0012305  0.0001221   10.08   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.374 on 389 degrees of freedom
## Multiple R-squared:  0.6876, Adjusted R-squared:  0.686 
## F-statistic:   428 on 2 and 389 DF,  p-value: < 2.2e-16

plot(predict(fit.lm), fit.lm$residuals, xlab = "Predicted", ylab = "Residuals")
abline(h = 0, col = 2, lty = 2)

• Remedy: 가중최소제곱법

3.7.3 Collinearity

• 다중공선성 문제 multicollinearity

  • 강한 상관관계가 있는 설명변수들을 모형에 포함시키는 경우 최소제곱추정에 문제를 발생시킴

  • 추정된 회귀계수값이 이상하거나 표준오차가 매우 커짐

  • 극단적인 경우 모형 적합이 아예 이루지지 않을 수 있음

pairs(Auto[,-(7:9)])

fit.lm <- lm(mpg ~ displacement + horsepower + weight + acceleration, data = Auto)
summary(fit.lm)

## 
## Call:
## lm(formula = mpg ~ displacement + horsepower + weight + acceleration, 
##     data = Auto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.378  -2.793  -0.333   2.193  16.256 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  45.2511397  2.4560447  18.424  < 2e-16 ***
## displacement -0.0060009  0.0067093  -0.894  0.37166    
## horsepower   -0.0436077  0.0165735  -2.631  0.00885 ** 
## weight       -0.0052805  0.0008109  -6.512  2.3e-10 ***
## acceleration -0.0231480  0.1256012  -0.184  0.85388    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.247 on 387 degrees of freedom
## Multiple R-squared:  0.707,  Adjusted R-squared:  0.704 
## F-statistic: 233.4 on 4 and 387 DF,  p-value: < 2.2e-16

• 분산팽창계수 VIF

  • Variance inflation factor (VIF) \[ \textsf{VIF}(\beta_j) = \frac{1}{1 - R_{j | -j}^2}\] 단, \(R_{j| -j}^2\)는 \(x_j\)를 반응변수로 하고 나머지 변수들을 설명변수로 해서 적합시킨 선형회귀모형의 \(R^2\)

  • 보통 VIF값이 10 이상인 설명변수가 있으면 다중공선성 문제가 있다고 봄

library(car)

## Loading required package: carData

vif(fit.lm)

## displacement   horsepower       weight acceleration 
##    10.686922     8.823022    10.284456     2.603255

• Remedies: 변수선택, regularizarion, PC regression, etc.

---

3.8 Logistic regression

3.8.1 Generalizations of the Linear Model

• Non-normal responses: Logistic regression, Poisson regression?

• Non-linearity: Kernel smoothing, splines and generalized additive models

• Interactions: Tree-based methods, bagging, random forests and boosting (these also capture non-linearities)

• Regularized fitting: Ridge regression and lasso

3.8.2 Logistic regression

• 지금까지의 회귀모형은 예측 대상 변수인 \(y\)값에 대해 정규분포 가정을 할 수 있는 경우였다. 그러나 실제 분석에서는 예측 대상 변수가 범주형이거나 계수형이어서 정규분포 가정이 적절치 않은 경우가 태반이다. 이러한 경우를 위한 선형모형을 일반화선형모형(generalized linear model, GLM)이라 한다.

• 로지스틱 회귀모형은 반응변수 \(y\)가 성공(1), 실패(0) 두 가지 값(예: 환자 생존 여부, 전염병 발병 여부, 신용 부도 발생 여부, 고객 claim 발생 여부 등)을 갖는 이항(binary) 변수인 경우에 대한 회귀모형으로 일반화선형모형의 대표적 방법론이다.

• 반응변수가 연속형인 경우와 달리, 이항 변수인 경우 \(E(y|x) = P(y=1\,|\,x)\)가 되어 취할 수 있는 값의 범위가 0과 1 사이의 값으로 제한된다. 이를 무시하고 통상적인 회귀분석을 적용하면 예측값이 허용 범위인 0과 1 사이를 벗어나는 경우가 발생해 문제가 된다.

• 따라서 다음과 같이 \({\rm logit} (p) = \log\frac{p}{1-p}\)을 이용(이런 함수를 link function이라 함)해 변환한 모형을 고려한다.
  \[ {\rm logit}P(y = 1\,|\,x_1,...,x_p) = \beta_0 + \beta_1 x_1 + ... + \beta_p x_p \] 참고로 로짓함수의 역함수는 sigmoid 함수 \(\sigma(s) = \frac{1}{1 + e^{-s}}\)이다. 즉, \[ P(y = 1\,|\,x_1,...,x_p) = \sigma(\beta_0 + \beta_1 x_1 + ... + \beta_p x_p) \] 으로 표현할 수 있다.

• 로짓함수 대신 표준정규분포 누적분포함수의 역함수를 사용하는 경우 프로빗(probit) 회귀모형이라 한다.

예제 데이터

• faraway 패키지의 built-in 데이터셋인 orings 데이터

  • 1986년 우주왕복선 챌린저호 폭발 사건이 로켓엔진의 O-ring 불량과 관련이 있다고 알려져있다.

  • 이 데이터셋은 이전 23회의 우주 왕복 임무 과정에서 수집된 데이터이다.

  • temp: 발사 때 기온 (화씨)

  • damage: 6회 시험 중 O-ring 손상이 발생한 횟수

library(faraway)

## 
## Attaching package: 'faraway'

## The following object is masked _by_ '.GlobalEnv':
## 
##     logit

## The following objects are masked from 'package:car':
## 
##     logit, vif

data(orings)
str(orings)

## 'data.frame':    23 obs. of  2 variables:
##  $ temp  : num  53 57 58 63 66 67 67 67 68 69 ...
##  $ damage: num  5 1 1 1 0 0 0 0 0 0 ...

fit.glm <- glm(cbind(damage, 6 - damage) ~ temp, data = orings,
               family = binomial(link = logit))
#fit.glm <- glm(cbind(damage, 6 - damage) ~ temp, data = orings, 
#               family = binomial(link = probit))
summary(fit.glm)

## 
## Call:
## glm(formula = cbind(damage, 6 - damage) ~ temp, family = binomial(link = logit), 
##     data = orings)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.9529  -0.7345  -0.4393  -0.2079   1.9565  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 11.66299    3.29626   3.538 0.000403 ***
## temp        -0.21623    0.05318  -4.066 4.78e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38.898  on 22  degrees of freedom
## Residual deviance: 16.912  on 21  degrees of freedom
## AIC: 33.675
## 
## Number of Fisher Scoring iterations: 6

plot(damage/6 ~ temp, xlim = c(25, 85), ylim = c(0, 1),
     col = "navy", pch = 20,
     xlab = "Temperature", ylab = "Pr(damage)", data = orings) 
x <- seq(from = 25, to = 85, by = 0.01) 
beta <- fit.glm$coefficients 
lines(x, sigmoid(beta[1] + beta[2]*x), col = 2, lwd = 3)
abline(h = 0.1, col = 2, lty = 2)
print(thre <- min(x[sigmoid(beta[1] + beta[2]*x) <= 0.1]))

## [1] 64.1

abline(v = thre, col = 2, lty = 2)
points(thre, 0, pch = 25, bg = "yellow", col = 2)

• 챌린저호 발사 당시 기온은 화씨 31도였다. O-ring 손상확률은 얼마인 지 예측해보자.

predict(fit.glm, newdata = data.frame(temp = 31), type = "response") # OMG!

##         1 
## 0.9930342